% Same as ../ppft/reconstruct3d2_interp, but uses the inversion of the PPFT
% that is based on a convolution filter.
%
% Yoel Shkolnisky, December 2010.

function vol=fastReconstruct3D(X,PHI,T1p,Tv,Nv,precision)

nt=size(X,1);
nr=size(X,2);
threadN=feature('numCores');
omega1_p=2*pi/T1p;
omega0_p=omega1_p/(2*nr+1);

T1v=Tv/Nv;
omega1_v=2*pi/T1v;

if (omega1_v<omega1_p)
    warning('Bandwidth of reconstructed volume is too small. Increase Nv');
end

pts=zeros(nr*nt,3,precision);
vals=zeros(nr*nt,1,precision);
for j=1:nr
    pts(j:nr:nr*nt,:)=PHI.*(j-1).*omega0_p;
    vals(j:nr:nr*nt)=X(:,j);
end

pp_pts1=zeros((Nv+1),(Nv+1),(3*Nv+1),3,precision);
pp_pts2=zeros((Nv+1),(Nv+1),(3*Nv+1),3,precision);
pp_pts3=zeros((Nv+1),(Nv+1),(3*Nv+1),3,precision);

pp_pts1_norm=zeros((Nv+1),(Nv+1),2,3,precision);
pp_pts2_norm=zeros((Nv+1),(Nv+1),2,3,precision);
pp_pts3_norm=zeros((Nv+1),(Nv+1),2,3,precision);

r=(-3*Nv/2:3*Nv/2).';

for l=-Nv/2:Nv/2
    for j=-Nv/2:Nv/2
        pp_pts1(j+Nv/2+1,l+Nv/2+1,:,:)=r*[1 (-2*l/Nv) (-2*j/Nv)]*omega1_v/(3*Nv+1);
        pp_pts2(j+Nv/2+1,l+Nv/2+1,:,:)=r*[(-2*l/Nv) 1 (-2*j/Nv)]*omega1_v/(3*Nv+1);
        pp_pts3(j+Nv/2+1,l+Nv/2+1,:,:)=r*[(-2*l/Nv) (-2*j/Nv) 1]*omega1_v/(3*Nv+1);
    end
end

pp_norm=zeros((Nv+1),(Nv+1),precision);

for l=-Nv/2:Nv/2
    for j=-Nv/2:Nv/2
        pp_norm(j+Nv/2+1,l+Nv/2+1)=sqrt(pp_pts1(j+Nv/2+1,l+Nv/2+1,1,1)^2+...
                            pp_pts1(j+Nv/2+1,l+Nv/2+1,1,2)^2+...
                            pp_pts1(j+Nv/2+1,l+Nv/2+1,1,3)^2);
    end
end

pp_pts1_norm(:,:,1,1)=pp_pts1(:,:,1,1)./pp_norm;
pp_pts1_norm(:,:,1,2)=pp_pts1(:,:,1,2)./pp_norm;
pp_pts1_norm(:,:,1,3)=pp_pts1(:,:,1,3)./pp_norm;
pp_pts2_norm(:,:,1,1)=pp_pts2(:,:,1,1)./pp_norm;
pp_pts2_norm(:,:,1,2)=pp_pts2(:,:,1,2)./pp_norm;
pp_pts2_norm(:,:,1,3)=pp_pts2(:,:,1,3)./pp_norm;
pp_pts3_norm(:,:,1,1)=pp_pts3(:,:,1,1)./pp_norm;
pp_pts3_norm(:,:,1,2)=pp_pts3(:,:,1,2)./pp_norm;
pp_pts3_norm(:,:,1,3)=pp_pts3(:,:,1,3)./pp_norm;

pp_pts1_norm(:,:,2,1)=pp_pts1(:,:,3*Nv+1,1)./pp_norm;
pp_pts1_norm(:,:,2,2)=pp_pts1(:,:,3*Nv+1,2)./pp_norm;
pp_pts1_norm(:,:,2,3)=pp_pts1(:,:,3*Nv+1,3)./pp_norm;
pp_pts2_norm(:,:,2,1)=pp_pts2(:,:,3*Nv+1,1)./pp_norm;
pp_pts2_norm(:,:,2,2)=pp_pts2(:,:,3*Nv+1,2)./pp_norm;
pp_pts2_norm(:,:,2,3)=pp_pts2(:,:,3*Nv+1,3)./pp_norm;
pp_pts3_norm(:,:,2,1)=pp_pts3(:,:,3*Nv+1,1)./pp_norm;
pp_pts3_norm(:,:,2,2)=pp_pts3(:,:,3*Nv+1,2)./pp_norm;
pp_pts3_norm(:,:,2,3)=pp_pts3(:,:,3*Nv+1,3)./pp_norm;

pp_pts1=reshape(pp_pts1,(3*Nv+1)*(Nv+1)*(Nv+1),3);
pp_pts2=reshape(pp_pts2,(3*Nv+1)*(Nv+1)*(Nv+1),3);
pp_pts3=reshape(pp_pts3,(3*Nv+1)*(Nv+1)*(Nv+1),3);

pp_pts1_norm=reshape(pp_pts1_norm,2*(Nv+1)*(Nv+1),3);
pp_pts2_norm=reshape(pp_pts2_norm,2*(Nv+1)*(Nv+1),3);
pp_pts3_norm=reshape(pp_pts3_norm,2*(Nv+1)*(Nv+1),3);

[az,elv,~] = cart2sph(PHI(:,1),PHI(:,2),PHI(:,3));

pp1_filter=zeros(nt,1,'int8');
pp2_filter=zeros(nt,1,'int8');
pp3_filter=zeros(nt,1,'int8');

for i=1:nt
    if ((abs(elv(i)) < pi/4+pi/32) && (abs(az(i)) < pi/4  +pi/32))
        pp1_filter(i)=1;
    end
    if ((abs(elv(i)) < pi/4+pi/32) && (abs(az(i)) > 3*pi/4-pi/32))
        pp1_filter(i)=-1;
    end
    if ((abs(elv(i)) < pi/4+pi/32) && (az(i) > pi/4-pi/32) && (az(i) < 3*pi/4+pi/32))
        pp2_filter(i)=1;
    end
    if ((abs(elv(i)) < pi/4+pi/32) && (az(i) < -pi/4+pi/32) && (az(i) > -3*pi/4-pi/32))
        pp2_filter(i)=-1;
    end
    if elv(i) > 7*asin(1/sqrt(3))/8
        pp3_filter(i)=1;
    end
    if elv(i) < -7*asin(1/sqrt(3))/8
        pp3_filter(i)=-1;
    end
end

knn=8;

if strcmp(precision,'single')
    [idx,dist1]=...
        KNN3DauxiliaryMEXSinglePrecision(threadN,nt,nr,Nv,1,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts1(:,1),pp_pts1(:,2),pp_pts1(:,3),pp1_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts1_norm(:,1),pp_pts1_norm(:,2),pp_pts1_norm(:,3));
    idx=idx+1;
    tmpvals1=interp(vals,idx,dist1);
    [idx,dist2]=...
        KNN3DauxiliaryMEXSinglePrecision(threadN,nt,nr,Nv,2,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts2(:,1),pp_pts2(:,2),pp_pts2(:,3),pp2_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts2_norm(:,1),pp_pts2_norm(:,2),pp_pts2_norm(:,3));
    idx=idx+1;
    tmpvals2=interp(vals,idx,dist2);
    [idx,dist3]=...
        KNN3DauxiliaryMEXSinglePrecision(threadN,nt,nr,Nv,3,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts3(:,1),pp_pts3(:,2),pp_pts3(:,3),pp3_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts3_norm(:,1),pp_pts3_norm(:,2),pp_pts3_norm(:,3));
    idx=idx+1;
    tmpvals3=interp(vals,idx,dist3);
elseif strcmp(precision,'double')
    [idx,dist1]=...
        KNN3DauxiliaryMEXDoublePrecision(threadN,nt,nr,Nv,1,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts1(:,1),pp_pts1(:,2),pp_pts1(:,3),pp1_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts1_norm(:,1),pp_pts1_norm(:,2),pp_pts1_norm(:,3));
    idx=idx+1;
    tmpvals1=interp(vals,idx,dist1);
    [idx,dist2]=...
        KNN3DauxiliaryMEXDoublePrecision(threadN,nt,nr,Nv,2,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts2(:,1),pp_pts2(:,2),pp_pts2(:,3),pp2_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts2_norm(:,1),pp_pts2_norm(:,2),pp_pts2_norm(:,3));
    idx=idx+1;
    tmpvals2=interp(vals,idx,dist2);
    [idx,dist3]=...
        KNN3DauxiliaryMEXDoublePrecision(threadN,nt,nr,Nv,3,knn,...
                        pts(:,1),pts(:,2),pts(:,3),...
                        pp_pts3(:,1),pp_pts3(:,2),pp_pts3(:,3),pp3_filter,...
                        PHI(:,1),PHI(:,2),PHI(:,3),...
                        pp_pts3_norm(:,1),pp_pts3_norm(:,2),pp_pts3_norm(:,3));
    idx=idx+1;
    tmpvals3=interp(vals,idx,dist3);
else
    fprintf('unknown precision %s\n',precision);
    vol = -1;
    return ;
end

ppvol=zeros(3,(3*Nv+1),(Nv+1),(Nv+1),precision);

idx=1;
for k=-3*Nv/2:3*Nv/2
    for l=-Nv/2:Nv/2
        for j=-Nv/2:Nv/2
            if dist1(1,idx) >= 0
                ppvol(1,k+3*Nv/2+1,l+Nv/2+1,j+Nv/2+1)=tmpvals1(idx);
            end
            if dist2(1,idx) >= 0
                ppvol(2,k+3*Nv/2+1,l+Nv/2+1,j+Nv/2+1)=tmpvals2(idx);
            end
            if dist3(1,idx) >= 0
                ppvol(3,k+3*Nv/2+1,l+Nv/2+1,j+Nv/2+1)=tmpvals3(idx);
            end
            idx=idx+1;
        end
    end
end

precond=1;

L=3*Nv-2;
filename=sprintf('filters%sH3Dfilter%dprecond%dprecision%s.dat',filesep,Nv,precond,precision);
if exist(filename,'file')
    H=load3DHFilter(filename,precision);
else
    H=save3DHFilter(filename,Nv,precond,precision);
end

Hhat=fftn(H.filter,[L L L]);

temp = preconditionedAdjointPPFT3D(ppvol,precision);
ii=norm(imag(temp(:)))/norm(temp(:));
fprintf('Imaginary temp part = %E\n',ii);

zerosLXLXL=zeros([L L L],precision);
zerosNvXNvXNv=zeros([Nv Nv Nv],precision);
[vol,~,~,~] = ConjugateGradient3D(temp,Hhat,Nv,zerosLXLXL,1.0e-10,60,zerosNvXNvXNv,1);

% If we take a volume, generate projections, and reconstruct, we want the
% reconstructed volume to be in the same spatial orientation as the
% original one. I discovered experimentally that this is achieved by
% flippining all dimensions. I haven't check why this happens, but in any
% case, we fix that by the following flips:

vol=flipdim(flipdim(flipdim(vol,1),2),3);

% This flip introduces a one pixel shift, so we fix that
vol2=zeros(size(vol),precision);
vol2(2:end,2:end,2:end)=vol(1:end-1,1:end-1,1:end-1);
vol=vol2;

end

function v=interp(vals,idxs,dists)
    mm=min(dists);  % Set the weight to the nearest neighbor to be 1/2.
    EPS=1.0E-12;
    mm(mm<EPS)=1; % Whenever the distance to the closest point is zero, set it arbitrarily to 1.
    delta=mm./log(1/2);
    w=exp(dists./repmat(delta,size(idxs,1),1));
    v=sum(vals(idxs).*w)./sum(w);
end
